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Distinct channels of interaction in a complex networked system define network layers, which co- 
exist and co-operate for the system's function. Towards realistic modeling and understanding such 
multiplex systems, we introduce and study a class of growing multiplex network models in which 
different network layers coevolve, and examine how the entangled growth of coevolving layers can 
shape the overall network structure. We show analytically and numerically that the coevolution 
can induce strong degree correlations across layers, as well as modulate degree distributions. We 
further show that such a coevolution-induced correlated multiplexity can alter the system's response 
to dynamical process, exemplified by the suppressed susceptibility to a threshold cascade process. 
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Introduction — Agents in complex systems interact in 
many ways: People are influenced by multiple chan- 
nels of social interaction such as friendship and work- 
partnership, and multiple means of transportation such 
as avian and ground transportations constitute the global 
transportation infrastructure [UJ3]. Such systems can 
best be represented as multiplex networks with multiple 
types of links. Each link- type in the system defines a net- 
work layer, which is by no means in isolation but co-exist 
and co-operate with other layers to fulfill the system's 
function. Such coupling and interplay of network layers 
can result in emergent structural and dynamical impact 
in nontrivial ways [4HTQ], rendering the understanding 
based on single- network approach incomplete. 

In real- world complex systems, either self-organized or 
man-made, the coupling between network layers is not 
completely random but structured, a property referred 
to as correlated multiplexity [5 . Specifically, degree of 
a node (degree is the number of links a node has) in 
one layer and that in the other are often strongly corre- 
lated, as evidenced in social pQ and transportation net- 
works [2] ITT] . Expectedly, she who has many friends 
tends to have many friendly coworkers in workplace; a 
hub-airport city is most likely a rail-hub, and so on. Such 
a structured coupling of network layers is shown to af- 
fect the system's connectivity and robustness properties 
[5] [TIHI3] . However, its underlying evolutionary mecha- 
nism has not yet been systematically investigated. 

In this paper we propose the coevolution of network 
layers as an evolutionary mechanism for the correlated 
multiplexity in growing multiplex networks. To substan- 
tiate the idea, let us turn back to the transportation net- 
work example. Suppose one were to establish a new air 
route. In doing so, it might be reasonable to consider not 
only the candidate city's avian connectivity, but also its 
ground connectivity such as rail and highway infrastruc- 
ture in order to maximize the synergy. That is, layers 
in a multiplex system do not merely co-exist, but they 
co-evolve, affecting each other's growth. Coevolution en- 
tangles the growth dynamics of different layers, and elu- 



cidating its effect on network structure and dynamics is 
the main aim of this paper. 

To this end, we introduce and study a class of grow- 
ing multiplex network models with coevolving layers. We 
show by analytic calculations assisted by extensive simu- 
lations that coevolution can profoundly affect the struc- 
ture of multiplex systems. Not only can it shape the cor- 
related multiplexity, it can also modulate the degree dis- 
tributions. We further demonstrate that coevolved mul- 
tiplex structure responds differently to a cascade process 
compared to independently evolved ones, exemplifying 
the dynamical signature that coevolution can imprint. 
Note that coevolution of (single) network structure and 
dynamical process on it has been studied [HUE]. Yet, 
coevolution effect of different layers within a multiplex 
system has remained unexplored. 

Coevolving multiplex network model. — To enlighten 
ourselves on the impact of coevolution, we build a sim- 
ple toy model of coevolving multiplex network (Fig. la). 
Each step, a new node enters into the system and in each 
layer establishes a link to an existing node. Probability 
that an existing node would receive a link from the new 
node gives the growth kernel II of its degree [16 . We 
formulate the coevolution of network layers in the way 
that the growth kernel of a node's degree in layer \i is 
not only dependent on its degree in that layer but also 
on its degrees in other layers, 

n M = f(k a , kp, . . . , k^, . . . , k e ), (1) 

where t is the total number of layers in the system and we 
used Greek subscripts to denote the layer index. We are 
interested in not only the degree distributions of layers 
grown under Eq. (1), but also the correlation of degrees 
across layers to address the correlated multiplexity in the 
multiplex system. For the latter, we calculate Pearson 
correlation coefficient p^ v between the degrees of a node 
in the two layers fi and v pQ , given as 

((fc-(fc))q-(Q)) = (ki)- (k)(i) 
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FIG. 1: (a) Model illustration. Each step, a new node enters 
the system and establishes a link in each layer. To choose the 
node to connect in the layer a, the new node refers to the net- 
work connectivity not only in that layer a but also in the other 
layer /3 (and similarly in the layer f3). Relative dependency to 
the other layer is controlled by the coevolution factor e. (b) 
Numerical simulation results for the joint degree distribution 
P(kct,kp) of coevolving networks of size N = 10 3 with sim- 
ple preferential attachment for various e. As the coevolution 
factor e increases, degrees of a node in the two layers become 
more strongly correlated, intensifying correlated multiplexity. 
Colorbar denotes the scale in \nP(k a , kp). 



where we used a simpler notation that k = k^ and I = k v 
and (k) (a^) is the mean (standard deviation) of node 
degrees in the given layer. 

Henceforth, we focus our analyses on a class of growth 
kernels based on linear preferential attachment [I6j [17] 
and systems with t = 2 layers, or duplex system, al- 
though the main messages would be applicable to more 
general cases. 

Mutually- dependent layers — Let us suppose the 
growth kernels for the two layers given by 

U a ex [(1 - e)(k a + a) + e(kp + a)], (3a) 
cx [e(k a + a) + (1 - e)(kp + a)]. (3b) 

Here e is a parameter that controls the strength of co- 
evolution, hence called the coevolution factor. If e = 0, 
two layers grow independently. If e > 0, they coe- 
volve, mutually depending more strongly as e increases. 
a is the shift factor introduced to control the layer's na- 
tive degree exponent. Recall that growing network with 
II(fc) cx (k + a) has an asymptotic power-law degree dis- 
tribution, P(k) ~ /c -7 , with the exponent 7 = 3 + a [IT] . 

The case with a = (simple preferential attach- 
ment [16 ) is particularly illustrative as it is amenable 
to most detailed analytic results. The rate equation for 



node z's degree in layer a take a simple form as 

dk^ot _ (1 - t)h,a + ekij 



dt 
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(4) 



and similarly for the layer (3. The solution of Eq. (4) takes 

the same form as the original Barabasi-Albert model as 

1/2 

ki(t) = (t/ti) for both layers, where U is the arrival 
time of node i [16]. This leads to scale-free network lay- 
ers with P{k) ~ k~ 3 for both layers, irrespective of the 
coevolution factor e. The degree correlation is, however, 
crucially affected by the coevolution (Fig. lb). 

To see this, we set up a rate equation for the number 
of nodes having degree k on the layer a and I on the layer 
f3 at time t, denoted as Ck,i(t), which reads 



dt 
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where the parenthesized superscript is used to denote the 
layer index. Changing variable by Ck,i(t) = tck,i(t) and 
introducing the generating function for Ck,i(t) [18 , one 
can obtain following coupled differential equations for 
(kl) and (fc 2 ). 
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e(kl) = e{k 2 )+t 
Solving for {k 2 ), one obtains 
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where Ei(x) is the exponential integral [19, and ci, C2 and 
d\ are constants determined by boundary conditions. 

When e = 0, Eqs. (6) decouple and from Eq. (6a) one 
obtains 



(kl) = exp 



1 



c 3 - Ei 



(8a) 



where cs is another constant determined by the boundary 
condition. For e > 0, by plugging Eq. (7) into Eq. (6b) 
we have 



(kl) = 21n(t) - d 



Ei 
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exp 
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(8b) 

Combining Eqs. (7) and (8), one can obtain from 
Eq. (2) the correlation coefficient p between degrees of 
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FIG. 2: Coevolution induces strong correlated multiplexity. 
(a) Degree correlation between two layers p as a function of 
coevolutin factor e and shift factor a. (b) Plots of ln(t) [1 — p(t)] 
as function of t with a = for different e. Horizontal behavior 
of the simulation results (points) supports the analytically- 
predicted logarithmic convergence of p towards the asymp- 
totic value 1. The height of horizontal lines gives the co- 
efficient c of logarithmic correction term, (c) Comparison 
between the numerical simulation for p (points) with the the- 
oretical result Eq. (10) (lines) as a function of e for various 
a, showing excellent agreements. Numerical simulations are 
performed with network size N = 10 6 , averaged over 10 3 runs. 



a node in the two layers. In the long time (equivalently, 
large network) limit, it is obtained as 



1/2 
1 



(6 = 0), 

(0 < e < 1). 



(9) 



Convergence to the asymptotic values are logarithmic, 
with the leading correction term ~ — c/\n(t) with c be- 
ing a constant dependent on e. Nonzero correlation even 
for independently evolving case can be attributed to the 
age effect inherent in the growing network [11]. Yet, with 
coevolution (e > 0), the correlation jumps to its maxi- 
mum possible value 1, yet with a logarithmically slow 
convergence (Fig. 2a). Such a logarithmic convergence 
to perfect correlated multiplexity is confirmed by numer- 
ical simulation (Fig. 2b). 

When a > (the shifted linear kernel), similar pro- 
cedure leads to P(k) ~ &-( 3 + a ) for both layers. Yet 
again, the coevolution factor can affect the correlation. 
As 7 > 3, (k 2 ) and (kl) remain finite as network grows 




FIG. 3: Coevolution can modulate the degree distribution. 
Plotted are the cumulative degree distribution Pp(> k) of the 
dependently-evolving layer f3 for different coevolution factor 
e = 0, 0.1, 0.9, 1.0 (top to bottom) with fixed a = 1. Straight 
lines indicate theoretical exponents from Eq. (12), with slopes 
—2, —20/9, and —3 (top to bottom). Data are obtained from 
the networks with size N = 10 7 , averaged over 10 2 runs. Note 
that the results for e = 0.9 and 1.0 almost overlap. 



and converge rapidly to the limiting value. It is thus suffi- 
cient to focus on the limiting values. Similar but slightly 
more involved calculations lead to 



6e + a 
6e + 2a' 



(10) 



in excellent agreement with numerical simulations 
(Fig. 2c). p increases with e. It decreases with a (or 
equivalently, 7), yet remains p > 1/2 as long as 7 is fi- 
nite (Fig. 2a). These results clearly highlight the role of 
coevolution factor in shaping correlated multiplexity. 

Unidirectional dependency and dissimilar kernel — 
Layers may influence each other asymmetrically and non- 
reciprocally. Furthermore, each layer may have different 
native growth dynamics. The former can be dictated by 
distinct e parameters for the two layers, and the latter by 
different a parameters. To illustrate the effects of such 
factors, we consider the growth kernel of the following 
form 



U a oc (k a + a), 

lip ex [ek a + (1 - e)kp] . 



(11a) 
(lib) 



That is, the layer a grows independently with shift fac- 
tor a, but the layer j3 evolves with coevolution factor e, 
representing cases with unidirectional dependency with 
dissimilar growth kernel. 

In this case even the degree equation becomes quite 
involved for layer /?, but the limiting behavior can be 
obtained that 



"(3-e)/(l- 
-(3+a) 



0) , 

1) , 



(12) 
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with the two regimes separated by e ~ a/ (2 -fa). For the 
independently evolving layer a, P a (k) ~ /c~( 3+a ). This 
result shows that the degree distribution of the dependent 
layer /3 becomes modulated by the degree distribution of 
the layer a it depends on, if the coe volution factor is 
strong enough. This intriguing analytical prediction is 
supported by numerical simulations (Fig. 3). 

Similarly, the asymptotic value of the degree correla- 
tion coefficient for coevolving case (e > 0) is obtained 
as 



(13) 



p again increases with the coevolution factor e. Note that 
when a = (unidirectional dependency with identical 
kernel), p approaches to 1 (again logarithmically slowly) 
as in the mutually-dependent case. 

Impact on cascade dynamics. — Finally, we study the 
effect of coevolved multiplex structure on dynamical pro- 
cesses occurring on it. As a specific example, we con- 
sider the binary-state linear threshold cascade model in- 
troduced by Watts [20] and recently generalized for mul- 
tiplex networks [7]. In this model, each node can be in 
either active or inactive state. In the original (single net- 
work) version, an inactive node switches to active state, 
if the fraction of active nodes among its neighbors ex- 
ceeds the prescribed threshold R. One is interested in 
the final fraction (j) of the active nodes in the network, 
starting from a small fraction </>o of initial active seed 
nodes, which measures how susceptible a network is to 
the cascade. In the multiplex version, a node gets acti- 
vated if the fraction of active nodes among its neighbors 
exceeds the threshold in any layer [7 . Here we take frac- 
tion 0o of highest degree nodes as initial active seeds, 
and measure what fraction <j> of nodes are activated at 
the end of the multiplex cascade process. 

To highlight the effect of coevolution, we first compare 
the cascade processes on independently evolved (e = 0) 
and strongly coevolved (e = 1) structures. Results for 
networks with layers growing under simple preferential 
attachment, Eq. (4), are shown in Fig. 4. For a wide 
range of threshold R, the coevolved multiplex structures 
support significantly smaller cascades than the indepen- 
dently evolved ones. For given i?, the cascade size mono- 
tonically decreases with the coevolution factor e (Fig. 4, 
inset), showing that the coevolved structure with strong 
correlated multiplexity can be significantly less suscepti- 
ble to cascades. Note that our result resonates with other 
recent studies on increased robustness of correlated inter- 
dependent networks under cascading failures [TTJ [12] and 
increased attack resilience of correlated multiplex net- 
works [13] . 

Summary — To summarize, we have studied the effect 
of coevolution of network layers in multiplex networks by 
introducing a class of coevolving multiplex network mod- 
els. We have shown both analytically and numerically 
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FIG. 4: Coevolution can suppress the network's susceptibil- 
ity to cascade process. (Main) Compared are the final cas- 
cade size as a function of the threshold R with </>o = 10 -2 , 
for strongly-coevolving (e = 1) (red solid) and independently 
evolving (e = 0) (blue dotted) networks of size N = 10 5 . (In- 
set) as a function of e for fixed R — 0.65, showing that (j) 
monotonically decreases with e. 



that the coevolution can profoundly alter the structural 
properties of the evolved network, both in the degree dis- 
tribution within the layer and in the degree correlation 
across the layers. Coevolved multiplex structures spon- 
taneously develop strong correlated multiplexity. Such 
a structural modulation of coevolved multiplex is fur- 
ther shown to entail dynamical signature, as exempli- 
fied by the suppressed susceptibility to a cascade pro- 
cess. Coevolution of layers or parts of the whole sys- 
tem takes place ubiquitously, from social systems pQ and 
transportation infrastructures [2] to economic [2Tj and 
ecosystems [22] . The current coevolution-based modeling 
framework will serve as a starting point for further inves- 
tigation in such diverse fields, on top of which system- 
specific contexts and details can be embedded. From a 
theoretical perspective, an important open question yet 
to be addressed is the role of antagonistic dependency 
between layers, an apt subject for future study. 
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